/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of UHM (Unassembled HyperMatrix) sparse direct solver.
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "uhm/common.hxx"
#include "uhm/const.hxx"

#include "uhm/matrix/uhm/matrix.hxx"
#include "uhm/matrix/uhm/fla.hxx"

namespace uhm {
  static inline int chol_piv_flat   ( int fs, int ss, 
                                      linal::Flat_ ATL, linal::Flat_ ATR,
                                      linal::Flat_ ABL, linal::Flat_ ABR,
                                      linal::Flat_ p );
  
  static inline int chol_piv_hier   ( int fs, int ss,
                                      linal::Hier_ ATL, linal::Hier_ ATR,
                                      linal::Hier_ ABL, linal::Hier_ ABR,
                                      linal::Hier_ p );
  
  static inline int chol_incpiv_hier( int fs, int ss,
                                      linal::Hier_ ATL, linal::Hier_ ATR,
                                      linal::Hier_ ABL, linal::Hier_ ABR,
                                      linal::Hier_ p );
  
  void Matrix_FLA_::chol_piv() {

#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix
    // ----------------------------------------------------------
    chol_piv_hier( this->fs, this->ss, 
                   this->hier.ATL, this->hier.ATR,
                   this->hier.ABL, this->hier.ABR,
                   this->hier.p );
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix
    // ----------------------------------------------------------
    chol_piv_flat( this->fs, this->ss, 
                   this->flat.ATL, this->flat.ATR,
                   this->flat.ABL, this->flat.ABR,
                   this->flat.p );
#endif
  }

  
  void Matrix_FLA_::chol_incpiv() {

#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix
    // ----------------------------------------------------------
    chol_incpiv_hier( this->fs, this->ss, 
                      this->hier.ATL, this->hier.ATR,
                      this->hier.ABL, this->hier.ABR,
                      this->hier.p );
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix
    // ----------------------------------------------------------
    fprintf(stderr, ">> Not support incremental chol for flat matrix\n");
    abort();
#endif
  }


  static inline int chol_piv_flat( int fs, int ss, 
                                   linal::Flat_ ATL, linal::Flat_ ATR,
                                   linal::Flat_ ABL, linal::Flat_ ABR,
                                   linal::Flat_ p ) {
    if (fs) 
      linal::host::internal::chol_piv( FLA_LOWER_TRIANGULAR, ~ATL, ~p );
  
    if (fs && ss) {
      linal::host::internal::apply_pivots( FLA_RIGHT, FLA_TRANSPOSE, ~p, ~ABL );

      FLA_Trsm( FLA_RIGHT, FLA_LOWER_TRIANGULAR, FLA_CONJ_TRANSPOSE,
		FLA_NONUNIT_DIAG, FLA_ONE, ~ATL, ~ABL );
      FLA_Herk( FLA_LOWER_TRIANGULAR, FLA_NO_TRANSPOSE, FLA_MINUS_ONE,
		~ABL, FLA_ONE, ~ABR );
    }
    return true;
  }

  static inline int chol_piv_hier( int fs, int ss, 
                                   linal::Hier_ ATL, linal::Hier_ ATR,
                                   linal::Hier_ ABL, linal::Hier_ ABR,
                                   linal::Hier_ p ) {
    if (fs) 
      linal::host::internal::chol_piv( FLA_LOWER_TRIANGULAR, 
                                       ~(ATL.flat()), ~(p.flat()) );
  
    if (fs && ss) {
      linal::apply_pivots( FLA_RIGHT, FLA_TRANSPOSE, p, ABL );

      linal::trsm( FLA_RIGHT, FLA_LOWER_TRIANGULAR, FLA_CONJ_TRANSPOSE,
                          FLA_NONUNIT_DIAG, FLA_ONE, ATL, ABL );
      linal::herk( FLA_LOWER_TRIANGULAR, FLA_NO_TRANSPOSE, FLA_MINUS_ONE,
                          ABL, FLA_ONE, ABR );
    }
    return true;
  }

  static inline int chol_incpiv_hier( int fs, int ss, 
                                      linal::Hier_ ATL, linal::Hier_ ATR,
                                      linal::Hier_ ABL, linal::Hier_ ABR,
                                      linal::Hier_ p ) {
    if (fs) 
      linal::chol_incpiv( FLA_LOWER_TRIANGULAR, ATL, p );

    if (fs && ss) {
      linal::apply_pivots( FLA_RIGHT, FLA_TRANSPOSE, p, ABL );

      linal::trsm( FLA_RIGHT, FLA_LOWER_TRIANGULAR, FLA_CONJ_TRANSPOSE,
                          FLA_NONUNIT_DIAG, FLA_ONE, ATL, ABL );
      linal::herk( FLA_LOWER_TRIANGULAR, FLA_NO_TRANSPOSE, FLA_MINUS_ONE,
                          ABL, FLA_ONE, ABR );
    }
    return true;
  }


}
